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Summary 


The scalar form of the approximate factorization method has been used to develop a new 
code for the solution of three-dimensional internal laminar and turbulent compressible flows. I he 
Navier-Stokes equations in their Reynolds-averaged form are iterated in time until a steady 
solution is reached. Evidence is given to the implicit and explicit artificial damping schemes that 
proved to be particularly efficient in speeding up convergence and enhancing the algorithm 
robustness. A conservative treatment of these terms at domain boundaries is proposed in order to 
avoid undesired mass and/or momentum artificial fluxes. Turbulence effects are accounted for b> 
the zero-equation Baldwin-Lomax turbulence model and the q-u> two-equation model. For the 
first, an investigation on the model behavior in the case of multiple boundaries is performed. The 
flow in a developing S-duct is then solved in the laminar regime at Reynolds number (Re) 790 
and in the turbulent regime at Re=40.000 using the Baldwin-Lomax model. The Stanitz elbow is 
then solved by using an inviscid version of the same code at Grid dependence and 

convergence rate are investigated showing that for this solver the implicit damping scheme may 
play a critical role for convergence characteristics. The same flow at Re=2.5-10 is solved with 
the Baldwin-Lomax and the q-w models. Both approaches show satisfactory agreement with 
experiments, although the q-u/ model is slightly more accurate. 


Chapter 1 
INTRODUCTION 

With the advent of more and more powerful supercomputers, the numerical solution of 
three-dimensional turbulent flows has become possible (ref. 1). Although it is well known that 
lower order turbulence closures fail to reproduce secondary motions of Prandtl’s second kind 
(turbulence driven), they normally succeed in predicting secondary flows of the first, kind (pressure 
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driven). Therefore, for many complex configurations it has been possible to obtain fairly accurate 
results with zero and two-equation turbulence models (ref. 2) with a reasonable prediction of 
pressure-driven secondary flows. 

For three-dimensional blade-passage flows, a correct prediction of the wake behavior has 
been obtained by Yokota (ref. 3) by means of the standard high-Reynolds-number form of the k-c 
two-equation model with the wall function approach. However the zero-equation model 
implemented in reference 3 did not give a satisfactory depiction of the flow in the wake region. 
For the prediction of blade pressure distribution, the Baldwin-Lomax (ref. 4) zero-equation 
turbulence model was found to give accurate results in several flow configurations (ref. 5) despite 
the low convergence rate. For incompressible internal flows with no separation, good results have 
been obtained by Towne (ref. 6) with a zero-equation turbulence model and a parabolized Navier- 
Stokes solver. Still, for practical flow configurations, it is necessary to face quite long computing 
times mainly because of the large number of points usually required for a detailed description of 
the flow field. Moreover, the non-linearities associated with nearly all the turbulence models can 
play a significant role in slowing down the convergence rate to the steady state solution. This 
behavior is independent of the implicit or explicit nature of the flow solver and is intrinsic to the 
turbulence models (zero- or two-equation.) 

In two dimensions, a wide variety of flow conditions have been accurately solved by means 
of low-Reynolds-number forms of the k-e model in which the effect of laminar viscosity is 
explicitly accounted for (see refs. 7-9). In nearly all the flow conditions investigated, these forms 
proved to be more accurate than the standard high-Re formulation, provided that a sufficient 
number of grid points are located inside the viscous and buffer layers; in fact, secondary mot ions 
and losses are mainly driven by what happens close to walls so that a correct description of this 
flow region is crucial for an accurate simulation of the flow pattern. Rodi (ref. 9) found that the 
use of the low-Reynolds-number forms of the k-t model could give the prediction of secondary 
flows normally lost with the high-Reynolds-number form. Unfortunately, the first author (ref. 7), 
found some of these forms extremely stiff from a numerical point of view. The stiffness was 
mainly caused by the low-Reynolds-number effect terms in which exponential functions are 
introduced to model the wall effects. From this standpoint, it appeared worthwhile investigating 
some features of the turbulence model for internal turbulent flows by using an implicit algorithm. 
Since complex flow patterns, such as separation and dominant viscous effects, are expected in 
internal flows, the implicit approach was selected to increase the robustness and convergence rate 
of the numerical procedure when zero- and two-equation turbulence models were used. 

Chapter 2 

DESCRIPTION OF THE ALGORITHM 

2.1 - Navier-Stokes Equations 

The Boussinesq hypothesis allows relating the turbulent shear stresses to the mean strains 
via the so-called “eddy viscosity" so that, under this assumption, the three-dimensional 
Reynolds-averaged Navier-Stokes (N-S) equations can be written in divergence form and, 
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subsequently, transformed from the Cartesian coordinate system (x,y,z) to the generalized 
curvilinear coordinate system £,q,C- The resulting set of equations can be written in vector form 
as follows: 


dQ , dE 
dt + di 


, gF , dG _ <9Ev , dF v dG v 
+ dr) + ” di + dr) + 


( 2 . 1 ) 


where the flux vectors are 
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U, V, W are the unsealed contravariant velocities defined as follows: 


t/= fx U + £yV + £ 2 W 
V — T7 X U + *7yV + 1) 2 W 
< X U + CyV + CzW 


The shear stresses in three dimensions are 



and 

1 ft 2 

E e = Ur xx -f Vr xy + ^ r xz + i le (y-T)~dx 

1 fta 1 2 

F e = Ur yx + Vr yy + Wr y2 -f 
G e = Ur zx + Vr 2y + W r 2Z + 

in which the following definition for the x, as well as for y and z Cartesian derivatives, holds. 


d __ c d 
fa “ ^*dl 


d_ 

K dv 


C — 


In this set of equations, p is the fluid density, p is the static pressure, e is the total energy 
per unit volume, a is the sound speed, and (U,V,W) are the Cartesian components of the velocity 
vector. According to the Boussincsq assumption, the diffusion coefficients for momentum and 
energy are defined as follows 


in which ^ is the laminar viscosity that was considered independent of the static temperature, 
and is the turbulent viscosity obtained from the turbulence model. In this set of calculations, 
the turbulent Prandtl number Pr t was set equal to 0.90 and the laminar Prandtl number Prj was 
0.72. The relation between static pressure and total energy is obtained through the equation of 
state for a perfect gas that yields 


P— (T'I) (e-p 


U 2 +V 2 +W 2 

2 


) 


in which 7=1.4 is the ratio of specific heat capacities for air. 


2.1.1 - Metric relations 

The metrics are usually obtained from a chain rule expansion of x^, x^, x^>, y^, y 7 p y^, z^, 
Zyj, z .. The following definitions of the metric relations hold: 

£x = j (y Tj z^-y^ Z77) ; Sy = J (z»; x^-z^ x^) ; £ y = J (y^ x T yi) x ; -) 

Vx = J (z^ y^-y^ z^) ; T 7 y = J (x^ z^-x^ z^) ; r/ z = J (y^ x^-x^ y^.) 

Cx = j (y^ Z^-Z^ y, ; ) ; <y = J (x jy z^-x^ z n ) ; < 2 - J (x^ y^-y^ x, ; ) 


and the Jacobian of the coordinate transformation is given by 

J' 1 = x^y^- + x^-y^z n + x^y^z^ - x^y^z n - x^y^ - x^y ^ 

The flux vectors are discretized by using centered finite differences. When the centered 
discretization is used for the metrics, it can be shown that in three dimensions the metric 
invariants are not satisfied (ref. 10). This may result in large discretization error. However, it is 
possible to satisfy the invariants by a simple averaging technique that gives metrics similar to 
those computed by a finite volume method. For example, £ x is computed as 

*$)->♦<(£)) 

in which \ is a central average operator. This averaging process was found to ensure better mass 
conservation properties in the present calculations, especially for highly stretched grids. 


2.1.2 - Nondimensionalization 

Since the solver is specifically designed to compute internal flows, the set of equations is 
nondimensionalized with respect to the total quantities in the inlet section, indicated with 
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subscript 01. 


U = = ; V = = ; W = w 

jaoi/r J a oi/T 

e - _ p - _ P . 

e ~ .2 / ’ P ~ Pol ’ P Poi ’ 

a 0l/7 



= 


£ 

^1 


Following this nondimensionalization, the flow Reynolds number (Re) is defined as 


Re = 


J(RToi) Poi^ 


in which L is typically a characteristic dimension of the inlet section. 


2.2 - Approximate Factorization; Scalar Form 

The approximate factorization method first proposed by Beam and Warming (ref. 11) splits 
an n-dimensional operator into the product of n one-dimensional operators. 1 his technique 
provides a strong link between the equations insofar as they are solved fully coupled. The main 
drawback to this method lies in the necessity of time-consuming block tn- or pentadiagona 
matrix inversion. This problem becomes more evident in three dimensions where the couplet 
solution yields a 5x5 block tridiagonal matrix. In order to make this algorithm more efficient and 
still maintain its strong implicit nature, Pulliam (refs. 10 and 12) proposed a scalar form of the 
approximate factorization. The form of the standard algorithm in three dimensions can be written 

as follows: 

[ I+eAt^A - <^A V )M I+QAt^TjB - ^B V )M I+eAt^C - <5^C V )]*AQ=RHS ^ 

in which I is the identity matrix, 0 is a parameter that allows weighting the explicit-implicit 
nature of the space operator (in the present calculations, since the steady state solution was 
sought, 0=1 was used that gave first-order accuracy in time), 6 is a centered difference operator, 
At is the time step, AQ=Q Yl+1 - Q n , and Q=J Q. The convective jacobians are A, B and C; A v , 
B v , and C v are the diffusive jacobians in the three directions £, r), (; and RHS represents the 
convective and diffusive fluxes defined as follows: 


A=m ; 

B=C ; 

c=f£ 

8Q ’ 

dQ 

<9Q 

A A ■ 
Av “3Q ’ 

r- 5F v 
1 B ~dQ 

p dG v 
; “ dQ 



RHS=At(- 


dE dF dG , dE v , <9F V , 5G V \ 
d( ' dv ' dC d£ + dr) + IK) 


The Jacobian matrices A, B, C, A v , B v , and C v , that can be found in Pulliam (ref. 10) are 
full and require a general inversion process. Pulliam and Chaussee (ref. 12) proposed a more 
efficient procedure by making some simplifications. First, the hyperbolic property of the 
convective jacobians allows the diagonalization as 


A = T ( T y 1 ; B=T f1 \ r) T^ 1 ; C = T f A c T^ 1 


(2.3) 


in which T represents the eigenvector matrices, defined in (ref. 12). The eigenvalues of the 
jacobians in three dimensions, A, are given by 


A f = D[U,U,U, (/+aJ^ 2 +^ y 2 +^ 2 ,t/-aJ^ 2 +^ y 2 +G 2 ] 

A r, = D[V,V, V, (2.4) 

A c = D[W, W, W, M/+aJc x 2 +C y 2 +Cz 2 ,^-aJ<: x 2 +Cy 2 +C 2 2 ] 


where D stands for the main diagonal only. It is now possible to introduce 2.3 into 2.2 with the 
eigenvalues defined in 2.4. By doing so, it is impossible to diagonalize both the convective and 
diffusive jacobians since they have a totally different set of eigenvectors. By neglecting the 
diffusive jacobians and assuming the eigenmatrix locally constant, equation 2.2 can be rewritten 
as 


^♦[I+eAt^A^A^I+eAt^A^/Ml+eAt^Ajjir^AQ^RHS (2.5) 

in which the two matrices, N=T^ Tjj and P=TJ^ T,, have the nice property of being solution 
independent so that they can be computed only once (ref. 12). Obviously, equation 2.5 is only an 
approximation of the full form but it requires only scalar tri- or pentadiagonal matrix inversion 
since the A matrices are diagonal; this implies a reduction of nearly 50 percent of the operations 
required by the standard algorithm. We note that, even though the interior domain is solved 
implicitly, the boundary conditions are imposed in a fully explicit manner. This was done setting 
AQ=:0 at the domain boundaries during the implicit sweeps. 


2.3 - Local Time Stepping 

Since the steady state solution is sought, the following local time-stepping formula based on 
the approximate constant CFL condition is introduced, in which the contribution of the diffusive 
terms is accounted for as follows: 
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At^ = |tf| + a Js 2 + 4y + £z + **eff Rel + $ + & 
At jj =11/1 + 8 Jrj| + rfi + f|i + ^ e ff Re' 1 (^x + 'Jy + »?i) 

At^ = |W| + a JCx + Cy + Cz + ^eff Rel (d + Cy + (*) 

At = CFL 

At^ + At,j + At^ 


2.4 - Implicit Treatment of Diffusive Terms 

For inviscid Hows, the only assumption made in deriving the algorithm is that the 
differentiation of the eigenmatrices is neglected in 2.5. Additionally, when so ving t e - 
equations, the jacobians of the diffusive terms are normally assumed to be negligible This does 
not cause any stability or convergence problems for external flows where the diffusion-dominated 
region is small and restricted to a limited part of the computational domain. But this 
simplification may cause troubles for internal flows where the diffusion-dominated region can be 
large. According to this, for internal flows it is convenient to introduce an approximation of the 
eigenvalues of the diffusive Jacobian and put it into equation 2.5. Two forms of this 
approximation have been considered. 


• Pulliam’s Approximation: 

Pulliam (ref. 10) proposes an approximate form of the diffusive jacobian eigenvalues that 
was obtained by intense numerical testing. This form is 


A^ v = (p p eff Re' 1 (£ x 2 +£y 2 +£z 2 ) ^ X ) ' 
A^ v = ( p Re’ 1 ( 0 x 2 + 0 y 2 + T ? 2 2 ) J X ) ' 
A^ v = {p p ef f Re" 1 (Cx 2 +Cy 2 +Cz 2 ) J *) * 


In the present set of calculations it was found convenient not to weight the eigenvalues with the 
Jacobian of the coordinate transformation J" 1 . These expressions are included on the implicit side 
of 2.5. For instance, the £ direction implicit operator reads 


[ I + 0 At («£ - 6^ Aj) ] 


• Present Approximation: 

The exact form of the diffusive jacobians can be computed from the related flux vectors. For 
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the three sweeps, the main diagonal of such a matrix may be conveniently approximated as 


A £ V = D [0> Q £, 7 Pr' 1 ^] 

A r/ = D [0,a^,o^,a^,7 Pr 1 a r ^] 

A v = D [0,a.,o.,a>,7 Pr _1 a J 

(2.7) 

in which 

<*{ = t*eft K * 1 J' 1 (^ 2 +e y 2 +ez 2 ) 

a r) = P e ff 1 J 1 ( , ?x 2 +Vy 2 + , ?z 2 ) ^ 

a c = /z eff Re ' 1 J ' 1 (G 2 +Cy 2 +Cz 2 ) 


Regarding the extradiagonal terms as negligible, the previous diagonal matrices are a good 
approximation of the diffusive term jacobians and can be put into the implicit side of 2.5. For 
instance, the £ direction implicit operator reads 


[I + © At ^ Aj) ] 

in which the diffusive terms contribution is approximated as a first-order derivative. 

It is possible to prove that the first approximation increases the main diagonal dominance 
by summing to it an extra term, while in the second approximation an artificial term is 
subtracted from the off diagonal components while leaving the main diagonal unchanged. A 
comparison of the two approaches, 2.6 and 2.7, was performed on a simple straight channel 
geometry at Re~1000, Mj n | et =0.3 in a laminar flow regime. For this very simple flow 
configuration there were no differences in convergence rate between the two approaches and it was 
also possible to drop the diffusive terms on the implicit side without altering convergence. 
Differences started to appear at Re=50 because of the highly diffusive nature of the flow. Figure 
1(a) shows the best convergence history of the algorithm without any implicit treatment of the 
diffusive terms that was obtained at CFL=5 (the lower curve refers to the averaged residuals, 
while the upper one to the maximum). Figures 1(b) and 1(c) refer to equations 2.6 and 2.7. There 
are no appreciable differences between the two convergence rates obtained at CFL=10 since both 
the curves show nearly the same slope. The same result could be obtained if the Baldwin-Lomax 
turbulence model were used. At least for this class of flows, the approximate implicit treatment of 
diffusive terms given by equation 2,7 did not prove to be more efficient than 2.6. Further testing 
is necessary in order to verify this result at higher Mach numbers when differences in the 
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convergence between 2.6 and 2.7 may appear since the density derivatives do not tend to vanish 
for compressible flows as they do for incompressible ones. 


2.5 - Addition of Explicit and Implicit Damping 

When using centered finite differences, it is necessary to introduce an artificial damping 
scheme in order to prevent odd-even velocity-pressure decoupling that occurs whenever the local 
Peclet number exceeds two in large gradients regions. Moreover, the three-dimensional 
approximate factorization method can be proved to be unconditionally unstable since, when 
performing a linear stability analysis, one of the amplification factors is found to be slightly larger 
than one. This weak instability can be easily overcome by using an artificial damping scheme. 
From the two-dimensional version of the approximate factorization, it is also known that fourth- 
order artificial damping is necessary to damp the numerical modes associated with the highest 
frequencies. Following the basic guideline given by Jameson et al. (ref. 13) and Pulliam (ref. 14) a 
blended implicit-explicit second- plus fourth-order nonlinear damping scheme was introduced in 
the present solver. 

• Explicit Second- and Fourth-Order Damping: 

In the original formulation, the artificial damping is simply equally scaled according to the 
local At in the three space directions. This could be done mainly because the damping scheme has 
been first applied to in viscid flows. For viscous flow calculations, Pulliam (ref. 10) found it 
convenient to scale the damping terms according to the directional spectral radius a . The scheme 
for the £ direction sweep yields to 


in w 


D^+D^ = <5^((<7J' 1 )j + ij i k+( aJ X )i j,k)( w U,k^£ ' w i j.k^(Qi j.k^ 

hicli the weights u/ 2 ^ and u/ 4 ^ are defined as 


( 2 . 8 ) 


T 


i.j.k 


|Pj+l,j,k ' 2 Pij.k + Pj-l,j,k 
|Pi+l,j,k + 2p i,j.k + Pi-1 j,k 


‘‘'hk — ^ ^ mal (T|-i j,k’ ^i,j,k’ ^i+lj,k) 


,( 2 ) 


k = *”<**(0., Cl''"*' At - Wy' k ) 


( 4 ) 


( 2 ) 


(2.9) 


and the unsealed spectral radius is defined as 


<Tj j k = |t/| + a ^ + £y + £z 
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Since the present set of calculations has been performed for shock-free flow fields, the shock 
sensor term defined in 2.9 was switched off, taking always Tj j k . Nevertheless, this modification 
does not really affect the damping scheme since the second-order term is active only in presence of 
shocks. Actually, this scheme gives only fourth-order damping in smooth shock-free flow fields so 
that in most of the calculations f2 was set equal to zero whereas the fourth-order weight, Q^ 4 , 
ranges from 1/16 to 1/128. Equation 2.8 is then added to the RHS of 2.5 

RHS = RHS + At (D^ 2) + D^ 4) + + D^ 4) + D^ 2) + D< 4) ) 


• Implicit Second- or Fourth-Order Damping: 

To enhance the algorithm stability and convergence rate, it is helpful to include the artificial 
damping terms on the implicit side. For both the second- and fourth-order damping, the implicit 
treatment augments the diagonal dominance of the scalar system with a beneficial effect on the 
convergence rate. As was pointed out by Pulliam (ref. 10), there is a stability limit for the 
weights that can be used for the artificial damping terms. This limit is actually connected to the 
magnitude of the amplification factor of the scheme modified with the artificial terms. To obtain 
the best convergence rates, the implicit damping had to be the exact jacobian of the explicit 
counterpart added to the right-hand side. If the second-order damping is treated implicitly, 
equation 2.5 must be modified to include the added implicit terms. For instance, the £ sweep 
implicit side will be 


{l+0At(4 A - «,(((.J-‘) 1+lJ , k +(<rJ-‘) u k ) (4IiMij. k )))} 

( 2 . 10 ) 

This form maintains the tridiagonal nature of the jacobian matrix. 

When the fourth-order damping is treated implicitly, the differentiation of the jacobian 
matrix gives a scalar pentadiagonal system that, for the £ sweep, can be written as 


{[+6At(< ? A { ‘ 

( 2 . 11 ) 

A brief set of tests on a straight channel geometry proved that the fourth-order implicit damping 
(eq. 2.11) could bring a large gain in convergence rate with respect to the second-order (eq. 2.10). 
For a straight three-dimensional channel with approximately 10 4 points with an inlet-section- 
width to channel- length ratio of 15, typical of internal flow geometries, the tridiagonal solver 
associated with equation 2.10 could be run at CFL=2, which gives the convergence history shown 
in figure 2(a), while the pentadiagonal solver associated with equation 2.11 was shown to be much 
more robust and gave the much higher convergence rate shown in figure 2(b) with CFL=10. 
Despite a 25 percent increase in computational time to invert the pentadiagonal matrix with 
respect to the three-diagonal one, the implict fourth-order option proved to be much faster and 
more robust, and was retained in all the calculations. 
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• Wall Treatment of Artificial Damping: 

It is well known that any kind of artificial damping term introduces an error that has to be 
minimized in order not to affect the solution heavily. Among the various approaches that can be 
found in the literature, the nonlinear damping formulation proposed by Jameson et al. (ref. 13) 
ensures that the second-order terms are introduced only at shocks while keeping the fourth-order 
in the smooth region of the computational domain. However, it is possible to prove that, if the 
presence of the boundaries is not properly accounted for when introducing the artificial terms, 
nonzero momentum and mass fluxes can be produced at the boundaries. This fact can be easily 
seen by considering figures 3(a) and (b). Figure 3(a) shows the fourth-order damping weights 
applied at every single point % for the £ direction. The fourth-order difference stencil used here is 

^(Qij.k) = Qi.2j.k-4 Qj_ 1J>k + 6 Q iJik -4 Qj +1J k + Q i+2J ,k 

The fourth-order damping is normally switched off at t=l and i=2 so that no special treatment at 
the wall is required. Figure 3(a) shows that in doing so the sum of the fourth-order damping 
weights for every location i is zero only in the internal flow field for *>5. This ensures that no net 
fluxes are added in the internal domain only, while the sum of the artificial damping weights 
yields to nonzero values for 1 <i< 4. These weights actually correspond to a nonzero third-order 
derivative centered at *=(2+1/2). The same procedure may be followed for the second-order 
damping that is switched off at i=l (fig. 3(b)). The sum of the weights for every location i is zero 
only for i> 3. This gives a first-order derivative centered at i=l/2 that corresponds to a first- 
order flux at the wall. These flux errors can be easily controlled in two dimensions by a grid 
refinement in the wall proximity where high gradients are expected. In three-dimensional internal 
flows we are forced to use coarser grids, and the wall boundary condition is applied on very large 
surfaces with the possible result of large mass errors. To control the mentioned flux errors, a 
third-order derivative, with the differencing stencil given by 

^(Qy.k) = Qi-l,j,k + 3 Qjj.k" 3 Qi+Ij lk + Qi+2 j,k 

was added at i=2 to balance the fourth-order damping weight. The procedure was followed to 
balance the second-order damping at the wall where a first-order derivative, with the differencing 
stencil given by 


— -Qi-X j.k + Qi,j,k 


was added at t=2. 
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Chapter 3 

TURBULENCE MODELS 


3.1 - Baldwin-Lomax for Multiple Boundaries 

The standard version of the Baldwin-Lomax (ref. 4) zero-equation turbulence model was 
implemented here. This two-layer model divides the flow field into an inner layer close to the wall 
(in which viscous effects are dominant) and an outer layer. The two layers are governed by the 
following set of equations and constants. 

• Inner Layer: 

inner = Re P M 
I = K y D 

D = (1 - exp (- y + /A + )) 
y + = y 4 Re w m ax 


in which Q is the vorticity, D is the Van Driest damping function, y is the wail distance and 
K=0.41 is the Von Karman constant. 

• Outer Layer: 

The outer viscosity is applied from the point at which /i t innef — /pouter 

^t, outer “ R e P Kc C cp F wake F k | efc) (y) 

in which K c =0.0168 is the Clauser constant and C cp = 1.6 is an empirical constant; and 

rwake =m ' n (y max ^ max ’ ymaxU d jf F max ) 

in which F max , U d|f and F k(eb are given by 

F(y) = y M D 

u dif = (-J U 2 +V 2 +W 2 ) max - (V+v 2 +w 2 ) min 

p — 1 

k,eb 1 4- 5.5 (C k | eP y / ymax) 

where C kjeb =0.3, and ymax is the wall distance at which F(y) is maximum. 
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This algebraic turbulence model was originally developed for a single boundary layer, but in 
three-dimensional internal flows the presence of multiple boundaries causes the interaction of 
more than one boundary layer. While the inner, or viscous, layer is driven by what happens on 
the closest wall only, for the outer layer an averaging procedure is necessary to account for the 
various wall effects. In the present set of calculations, three approaches have been examined to 
account for the presence of more than one boundary layer. 

• Wall Treatment 1: 

Only the geometrically closest wall is considered for computing the outer viscosity without 
any averaging. 

• Wall Treatment 2: 

The inner layers are driven by the closest wall only, while the viscosity in the outer layer is 
computed with a simple weighted average according to the inverse of each wall distance as follows 


_ 1 ^t.outer^ 

^t. outer - n^y- L, l W w ’ 

VWw 

where W w is the wall distance, w is the wall number, and N is number of walls present in a cross 
section. 

• Wall Treatment 3: 

If only one boundary layer is present, the Van Driest damping succeeds in modeling the wall 
effect. Starting from this standpoint, the inner layer viscosities are computed by using only the 
closest wall contribution, while the viscosity in the outer layer is computed as a simple weighted 
average according to the inverse of the value of the Van Driest damping expression for each wall. 


3.2 - q-w Two-Equation Model 

In a previous investigation, Michelassi (ref. 7) found that the low Reynolds number forms of 
the two-equation models, such as k-€, could give an accurate prediction of two-dimensional 
incompressible separated flows. Unfortunately, these forms were found to be numerically stiff, 
mainly because of the correction terms introduced to model the low Reynolds effects that 
necessitate a strong mesh refinement in the sublayer. Furthermore, an initial profile for the 
turbulent quantities must be specified consistently to start the calculations. A first attempt to 
implement the Chien and the Rodi two-layer low Reynolds number forms of the k-6 model did 
not bring any converged result mainly because of difficulties in specifying both the mesh 
refinement and initial profiles for complex three-dimensional flows. Coakley (ref. 15), 
reassembling the Jones and Launder low Reynolds number form of the k-e model, proposed the q- 
u) two-equation model, in which the effect of molecular viscosity is directly modeled. This 
formulation ensured a better numerical behavior as compared with other low Reynolds number 
formulations. This model, rewritten in our notation, is 
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in which 




D = 1 - e' 


a R 


R = Re p q y 


p d = 


dUj 

3xj 




•I 4 u 


p 


\d _ Uj 
d / 5xj 


P = fi t S -l k P d 
C x = 0.405 D + 0.045 

= 0.09 ; C 2 = 0.92 ; a = 0.0065 ; <r q = 1.0 ; c t = 1.3 


The turbulent viscosity is computed as 

D ~~j- 

The two transported quantities, q and u>, are related to the more familiar k and e via the 
following relations: 

q = k 1/2 ; w = e / k (3-2) 


It is important to observe that e, appearing in equation 3.2, is the isotropic part of the 
dissipation rate. This quantity does not account for any nonisotropic effect (for example, the 
presence of a wall) and tends to zero on solid boundaries. (Conversely, the total dissipation rate 
tends to a finite value related to the wall shear stress.) This choice allows use of cj as an unknown 
since, by assuming that both k and i are going to zero at the wall with the same slope, tends to 
a finite value. 


3.2.1 - q-qj Solution 

The two transport equations for q and cj are implicitly solved with the same algorithm 
given in equation 2.2. The two equations are solved in a sequential manner and decoupled from 
the flow variables mainly because the coupling is provided only by the diffusive terms coefficients 
and the sink and source terms. Due to this choice, the scalar tridiagonal algorithm was 
implemented for the turbulence model solution. The only difference with respect to the solution of 
the N-S equations is the presence of the sink-source vector H. This term can be included in the 
three sweeps: 
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[ I+0At(-Hj c^+ 6 ( A - 6^A V )]* 

[ I+0At^-Hj Cjj + SjjB - 
[ I+0At(-Hj C( , + S^C - ^C V )]*AQ = RHS 


(3.3) 


where the same definitions given for equation 2.2 hold with only the addition of the jacobian of 
the sink and source terms, Hj, that is weighted in three sweeps according to c^, c^, c^ 
(c^+ Ct? +c£= 1). This jacobian is computed by neglecting the contribution of the damping 
function U, and for the two q and u equations reads 


h j 


d H q _ 1 DS 1 D u 
d(py\) "2 “ - 3 - 2 


w _ d H 1 ^ _ 2 r 

j d(pY l w) "3 1 


D P d - C 2 


2w 


In place of their exact form, Coakley proposes an approximation of the jacobians based on 
the turbulent viscosities that should ensure the dominance of the main diagonal. Figure 4 shows 
the comparison of the convergence rates of every single variable obtained without any sink or 
source terms jacobian, with the exact jacobian, and with Coakley’s approximation in a typical 
internal flow geometry. Surprisingly, there is no big gain in introducing the jacobian in the 
implicit side of the operator. The choice of the sweeps in which the two jacobians and Hj^ are 
introduced is not important. The error introduced in the approximate factorization of the implicit 
side of 3.3 increases roughly by a factor proportional to H: when introducing the jacobian in the 
three sweeps, thereby choosing C£=c^=c^ = l/3. This has only a weak influence on the 
convergence rate. Nevertheless, in tne present calculations the best convergence rates have been 
typically obtained by using c^=0, c ? ^=0.5, and c^=0.5, where £ is the main flow direction and 
T]X are the fine grid directions. 

While physical evidence shows that the turbulent kinetic energy k is zero at solid walls, the 
boundary condition for e, and consequently u>, is less evident. For the q-w model, Coakley (ref. 
15) found it convenient to impose a zero-normal derivative at the wall; this condition was 
retained in the present calculations. 
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Chapter 4 
RESULTS 


4.1 - Incompressible S-duct 

A first validation of the code was performed by computing the flow in an S-duct with a 
constant area cross section. The measurements (ref. 16) were taken for an incompressible fluid 
(water). The S-duct was given by two 22.5 degree bends with a 4-cm hydraulic diameter and 28- 
cm mean radius of curvature. This geometry was regarded as an interesting test since flow 
passages with similar shapes are often used to redirect the flow for air intakes in aeronautical 
engines. The efficiency of such ducts may be heavily affected by the presence of secondary flows so 
that the ability of a code to detect secondary velocities is essential for a proper depiction of the 
flow pattern. 

A sketch of the flow domain is shown in figure 5. Since the geometry and the flow were 
symmetric with respect to the x-r plane, it was possible to study only one-half of the duct 
imposing a symmetry boundary condition on the same plane. For the laminar flow regime, the 
Reynolds number, Re, based on the inlet bulk velocity, Ub, and hydraulic radius was 790. The 
grids employed for this calculation are shown in figures 6(a) and (b). The coarser grid (fig. 6(a)) 
has 70x29x15 points with a ratio between two consecutive grid cells in the cross-flow directions 
of 1.1 while the refined one (figure 6(b)) has 80x59x30 points with the same stretching ratio in 
the cross-flow plane. The fluid adopted in the experiments was water so that, in order to have 
negligible compressibility effects, the isentropic Mach number was set equal to 0.1 (that gives the 
inlet total pressure, outlet static pressure ratio) to ensure minor density changes. 

The flow pattern and the growth of the secondary motions are mainly pressure-driven 
because of the very smooth-bending walls that cause no flow separation. Still, this kind of flow 
necessitates a very accurate prediction of the boundary layer, even in the laminar regime, 
otherwise the secondary flows may be incorrectly predicted or completely lost. Figure 7(a) shows 
the measured (circles) and computed velocity profiles with the two grids mentioned above for the 
five sections, indicated in figure 5, on the symmetry plane reported in the experiments. The slight 
discrepancies in the computed velocities with the two grids at section 1, where the double-S starts, 
may be attributed to the thinner boundary layer predicted with the refined grid that produced a 
flatter velocity profile. Nevertheless, the agreement with experiments is fairly good. The 
agreement does not deteriorate for sections 2, 3, and 4 in either the coarse or refined grids. Section 
5 clearly shows that while the secondary motions could not be predicted by the coarse grid, they 
are fairly well reproduced by the refined one. This behavior is thought to be independent of the 
grid points in the main flow direction (as will be demonstrated in a further test) in which only 10 
points are added in the refined grid, and is closely related to the poor resolution of the boundary 
layer provided by the coarse grid in which the cross-flow momentum diffusion is clearly 
overestimated. 

Figure 7(b) shows the velocity profiles for the midspan section (r=l/2) at the same five 
sections. Basically the same comments made for the symmetry plane could be repeated here for 
the first three sections, while on sections 4 and 5, agreement deteriorates in the proximity of the 
symmetry plane. This may l>e attributed not only to the poor grid quality close to the symmetry 
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plane, but also to the fact that in the calculations a zero gradient condition was imposed in the 
direction normal to that plane, while experiments show that for section 4 the velocity gradient is 
far from being zero on the same plane. The typical computed secondary motion pattern at the 
exit of the second bend is shown in figure 7(c). The complex flow pattern exhibits two 
counterrotating vortices in the proximity of the two curved walls. 

Nearly 1500 iterations were necessary to reach an averaged residual of the order of 10 5 on 
the refined grid, while typically less than half that number of iterations were necessary to obtain 
the same residual with the coarser grid. These calculations were performed with an early version 
of the code where only the second-order implicit damping was implemented, so that the slow 
convergence rate may be attributed to both the small Mach number and to the small CFL 
number that could not exceed 2 without encountering stability problems. 

The turbulent flow regime was run at the experimental condition of Re=40.000. For this 
calculation, it was necessary to provide a more stretched grid at the wall, so that the 80x59x30 
grid was reassembled with a point expansion ratio in the cross-flow directions equal to 1.3. This 
provided the necessary point clustering at the boundaries to describe the thin boundary layer and 
allowed placing the first point at the wall at y^"ss2. For this test case, the Baldwin-Lomax model 
was implemented with wall treatment number 1. (See the section on Baldwin-Lomax for multiple 
boundaries.) 

Figure 8(a) shows the set of measured and computed velocity profiles on the symmetry 
plane for the same five cross sections given in figure 5. For this flow configuration, the agreement 
is reasonably good, especially at cross sections 4 and 5. In these sections the pressure gradients 
induce a strong secondary motion that seems to be correctly predicted by the present solver 
insofar as the agreement with measurements does not deteriorate as the second bend exit is 
approached and the secondary velocities reach their maximum. The momentum transfer to the 
external part of the second bend, particularly evident in sections 4 end 5, is well reproduced since 
the computed velocity profile asymmetry seems to be in close agreement with experiments. 

Figure 8(b) shows the computed and measured velocities on the midspan plane along the 
duct. The agreement is again good for the five sections and better than that found for the laminar 
flow regime mainly because the measured velocity profiles exhibit a zero-gradient on the 
symmetry plane that is correctly modeled by the symmetry condition imposed on the x-r plane. 
Still, at sections 3 and 4, the kink in the velocity profile close to the wall is not correctly 
predicted. 

The static pressure coefficient in the main flow direction, defined as 


C p - 


P " Pjniet 


is shown figure in 8(c) at three different locations (z=0,r=0), (z=0,r=l), and (r— 1/2, z— 1/2). Here 
the agreement is generally good. The pressure trend is correctly predicted together with the head 
loss. 


Figure 8(d) shows the convergence history obtained with CFL=2. For this calculation, the 
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same early version of the code mentioned above was used. The solid upper line refers to the 
maximum residual while the dashed line refers to the averaged residual; the spikes present in the 
first curve correspond to the updating of the turbulent viscosity performed every five iterations. It 
is remarkable that the pressure distribution did not change after the first 500 iterations, while to 
get the correct velocity profiles, it was found necessary to reach a residual of the order of 10 to 
10 ~ 6 . 


4.2 - Stanitz Elbow 

The flow in the accelerating rectangular elbow with 90 degree turning and the variable cross 
section described by Stanitz (ref. 17) has been computed with an inviscid version of the code and 
with both the Baldwin-Lomax and q -u turbulence models. This test case was selected because it 
provides a good set of measurements, including wall pressure distribution and visualization of 
secondary flows at the elbow exit section. The shape of the elbow was analytically computed by 
Stanitz (ref. 17) to give no separation with a strong area reduction and a specified pressure 
distribution on the side wall under incompressible flow conditions. A sketch of the experimental 
setup is shown in figure 9. Both the flow and the elbow geometry are symmetric with respect to 
the x-y midspan plane, thereby allowing a zero normal gradient condition. Among the various 
flow conditions investigated in reference 17, the one with M exjt =0.4 and with no spoiler at the 
duct inlet with a thin initial boundary layer was selected. 


4.2.1 - Inviscid Calculations 

The first set of tests was performed with an inviscid version of the code. T he convergence 
characteristics of the scalar form of the approximate factorization could be tested for inviscid 
calculations where the necessity of accounting for the diffusive terms on the implicit side of 
equation 2.5 drops. Only the pressure distribution on the walls of the elbow was compared with 
measurements in this set of calculations; no attempt was made to specify an experimental inlet 
profile of total pressure that was kept flat. The inlet boundary condition was specified by 
extrapolating the Riemann invariant from the first section inside the duct (ref. 18). Two grids 
were used: a 25x15x11 and a 50x15x11 with constant grid spacing in the cross-flow directions. 
The grid point locations in the streamwise direction were made to coincide with the points 
supplied by Stanitz (ref. 17) for the description of the elbow geometry, with the addition of three 
sections at the outlet to allow using a zero gradient condition. 

Figures 10(a) and (b) show the qualitative static pressure and Mach number isolines on the 
symmetry plane of the elbow. The small wiggles visible at the domain exit are due t ( o } the very 
small fourth-order damping weight that was set equal to ft =1/256, together with ft -0. The 
small value of the fourth-order damping weight is allowed by the very coarse grid in the cross- 
flow direction that automatically introduced a numerical diffusion. Despite the very coarse grids 
implemented here, the static pressure distribution p$, defined as 


P _ P ~ Pexit 
s Ptotal ' Pexit 
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and shown in figures 11(a) and (b), reveals fairly good agreement with experiments. While the 
pressure drop position is located correctly for the two grids on both the side and symmetry walls, 
the kink in the static pressure distribution on the suction surface at the side wall could not be 
reproduced since it was caused by the presence of secondary flows. The pressure distribution 
appears to be independent of secondary flows induced by viscous effects until the first part of the 
bend is reached. The local pressure rise located at S=2, where S is the streamwise coordinate 
along the centerline, is introduced by the growth of secondary velocities and is totally lost by the 
inviscid solver. The solution proved also to be fairly grid-independent, at least to the grid 
refinement in the main flow direction. No additional tests were performed to verify the influence 
of the point distribution in the cross-stream direction. 

The implicit treatment of second- and fourth-order artificial damping was compared by 
using the 50x15x11 grid. Figure 12 shows the comparison between the two implicit damping 
schemes. The solid line refers to the second-order implicit damping with ft (2 =1/4, while the 
dashed line refers to the fourth-order implicit option with ft =1/256. The gain in convergence 
rate is remarkable, in fact, the fourth-order solver could be run at CFL=10 while the second-order 
one could not be run at CFL>5. Any futher increase of the artificial damping weights gave slower 
convergence histories. Unfortunately, the same convergence rate is not obtainable for viscous 
calculations because of the strong point clustering and viscous effects that are not exactly 
accounted for on the implicit side of the operator. 


4.2.2 - Viscous Calculations 

The viscous calculations in the turbulent flow regime were performed on the five grids 
summarized in table I. The use of various point clusterings allowed a comprehensive investigation 
into the mesh dependence of the calculations. With this set of grids, it was possible to verify the 
influence of the grid point numbers in the main flow direction (with 51 or 99 points) and the 
cross-flow direction (with 31x21 and 41x31 points) with expansion ratios of 1.2 and 1.3. The 
refined grid (number 5) shown in figure 13 adopts the same point distributions in the main flow 
direction that were used for the refined grid in the inviscid calculations, and allows placing the 
first grid point at the wall at y ^1. The Reynolds number, based on the total conditions at the 
inlet section, is approximately 2.5- 10^~ 6 . 

These calculations were mainly aimed at the proper prediction of the wall pressure 
distribution that is heavily affected by the growth of secondary velocities. The choice of the 
experimental spoilerless configuration allowed comparing the turbulence models for a very thin 
boundary layer that required a heavy point stretching at the wall. Regarding the inlet boundary 
condition, in the Baldwin-Lomax model the inlet turbulent viscosity was extrapolated from inside 
the domain, and in the q-u; model a flat turbulent kinetic energy profile with various turbulence 
levels was specified at the inlet section, while w was extrapolated from inside the domain. 

The first set of tests concerned the Baldwin-Lomax model with the different wall treatments 
mentioned in section 3.1. With the experimental total pressure profile specified at the inlet 
section, the computed static pressure profiles are compared with the measurements in figure 14. 
Ihe plots refer to the static pressure distribution in the section corners on the side wall and the 
symmetry plane. It is evident that the way the outer viscosity is computed may play a significant 
role in the correct prediction of the pressure distribution. Wall treatment number 1, in which only 
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the closest wall is considered, and number 2, in which the outer turbulent viscosity is weighted 
according to the inverse of the wall distance, do not show large changes even if the two 
approaches are considerably different. For both techniques, the agreement with experimental 
results is fairly good on the pressure side of both the side wall and the symmetry plane. The 
suction side on both planes shows that the static pressure is overestimated. This is probably due 
to the presence of computed secondary flows that are much stronger than the experimental ones. 
The computed pressure drop induced by the presence of the bend and the strong flow acceleration 
is smoother than the measured one. This phenomenon is caused by high turbulent viscosities that 
induce a heavy momentum diffusion in the cross-flow direction, followed by a static pressure 
redistribution and growth of the boundary layer thickness. The agreement with measurements 
improves with multiple wall treatment number 3. No great changes between the multiple wall 
treatments are found for the pressure side of both the side wall and the symmetry plane where it 
was always possible to have accurate results. Still, the pressure minimum located at S=2, where 
velocities on the cross section start to develop, is better reproduced by wall treatment number 3. 
This test suggests that the averaging technique based on the Van Driest damping expression for 
the mixing length can give reasonable predictions. All the computations with the Baldwin-Lomax 
model were performed with multiple wall treatment number 3. 

The results of the mesh dependence tests for the Baldwin-Lomax model are summarized in 
figure 15 where the computed static pressure distribution in the four corners of the cross sections 
using different grids are compared with measurements. The pressure distribution profiles show 
that there are practically no differences between the predictions obtained with grids number 1 and 
3. This proves that an increment of the grid point numbers in the main flow direction does not 
produce any gain in terms of accuracy. This is strictly connected to the boundary layer resolution 
that is not improved by using grid number 3. The secondary flow growth is influenced by the low 
momentum regions located close to the wall, the correct simulation of which is not ensured by the 
two grids. Conversely, the implementation of grid number 2 clearly improves the accuracy of the 
results. The static pressure profile on the suction side of the side wall is in better agreement with 
experiments than the pressure distributions given by grids number 1 and 3. It is worthwhile 
observing that the use of a more refined grid in the cross-flow direction shows that the computed 
pressure minimum on the suction side is correctly located, even if its value is still overestimated, 
while this local minimum, located approximately at S=2, is completely lost with the other two 
grids. Nevertheless, the static pressure distribution on the symmetry plane appears to be weakly 
affected by grid refinement. No further investigation was performed by varying the cross-flow 
points expansion ratio. 

The flow simulation with the q-u> model required more tests since it has been necessary to 
investigate the dependence on both the mesh refinement and on the inlet turbulence level. The use 
of grids with an expansion ratio equal to 1.2 did not ensure significant improvements in results 
since not enough points were located close to the wall. In order to have a reasonable definition of 
the turbulent kinetic energy peak at the wall, the grid expansion ratio was fixed at 1.3. The 
results of the first set of tests are summarized in figure 16 in which, by adopting the coarse grid 
number 4, the inlet turbulence level was changed to verify its influence on the static pressure 
distribution. When decreasing the inlet turbulence level from 5.0 to 0.1 percent, the computed 
profiles progressively approach the measurements. Still, the predictions are far from the 
experiments for the suction side of both the symmetry plane and the side wall. This indicates the 
presence of a large momentum diffusion that is possibly caused by insufficient mesh refinement or 
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too high a turbulence level, or both. Figure 17 shows the results obtained with grid number 5 and 
lower turbulence levels. Figures 16 and 17 show that the refined grid with the same turbulence 
level brings some improvement in the agreement with the experiments, proving that the coarser 
grid was largely inadequate for this flow configuration (see, for instance, the 0.1 percent level). 
The growth of the secondary flows is evident in figure 17 where the static pressure distribution 
computed with the inviscid approach is compared with the results of the two-equation model 
obtained with different turbulence levels. 

A direct comparison of the turbulent calculations with the inviscid computation, in which 
the same pressure distribution is found on both the side wall and the symmetry plane because of 
the absence of secondary motions, shows large differences on the suction side only. These 
differences start from S=2 where experiments deviate from the inviscid solution. The one percent 
turbulence level, not reported in figure 17, appeared to bring too high a momentum diffusion. 
Consequently, this level was progressively decreased from 0.5 to 0.1 percent. The static pressure 
distribution on the pressure wall is mainly driven by convective phenomena, while the 
distribution on the suction wall is largely influenced by diffusion processes; this is why the 
pressure side distribution is always well reproduced and is nearly independent of the inlet 
turbulence level. The final result obtained with grid number 5 and a 0.1 percent turbulence level 
shows a fairly good agreement with experiments on both the suction and pressure sides, and seems 
to reproduce quite correctly the location and influence of secondary flows. 

The predicted pressure distributions obtained by the Baldwin-Lomax model with wall 
treatment number 3 and the q-u; model were very similar, but the two-equation model proved to 
be marginally more accurate, especially because of the static pressure distribution on the suction 
side. Figure 18 shows that the two models predict approximately the same static pressure pattern 
on the symmetry plane with nearly the same pressure drop due to the acceleration on the flow. 
Still, from figure 19, where the Mach number isolines are shown on the same symmetry plane, it 
is possible to observe that the zero equation model predicts a slightly thicker boundary layer than 
the one predicted with the q-u;. These differences start immediately after the inlet section and 
become more evident as the exit section is approached. This indicates a different turbulent 
viscosity distribution in the wall region. The Baldwin-Lomax model was in fact found to predict a 
sharper growth of the turbulent viscosity in the wall region. Still, the two models gave 
comparable values of the turbulent viscosity in the flow core. The differences in the boundary 
layer thickness are evident in figure 20 where the velocities in the exit section of the channel are 
plotted. While both turbulence models show the same flow pattern, the two-equation model 
predicts the center of the secondary recirculation closer to the wall than that predicted by the 
zero-equation model, in which the location of the center appeared to be in better agreement with 
experiments. This confirms a weaker cross-flow momentum diffusion (given by the two-equation 
model) as compared with the zero-equation. This is caused by the aforementioned differences in 
the turbulent viscosities. It is remarkable that both formulations predict a small secondary vortex 
in the wall corner of the suction and pressure walls. 

An interesting qualitative comparison of the predicted and measured secondary flows is 
given by figure 21. Figure 21(a) shows the experimental flow visualization by injecting a smoke 
filament inside the boundary layer close to the side wall at the inlet section of a reduced model of 
the channel. The low-momentum particles located well within the boundary layer are bent toward 
the suction wall by the pressure gradients, while the high-momentum particles exhibit a weaker 
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turning. Figure 21(b) shows the secondary flow predicted by the zero equation model, while figure 
21(c) refers to the results obtained with the q-w model. From this comparison, it is evident that 
the zero-equation model predicts a smoother turning of the velocities. This is in slightly better 
agreement with the experimental picture than the q-w model. 

A typical total pressure loss distribution on the channel exit section is given in figure 22. 
The results computed with the Baldwin-Lomax model using grids number 1 and 2 show the much 
thicker boundary layer obtained with the coarser grid. The refined grid ensures a correct 
description of secondary flows together with the location and magnitude of the total pressure 
losses. 

With the present research version of the code, 6000 sec were necessary to perform 2500 
iterations with the q-u> turbulence model with a 63.550 point grid on a CRAY-YMP 
supercomputer. This ensured an overall residual of the order of 10’ 6 , while a reduction of 
approximately 30 percent in CPU time could be obtained by using the zero-equation model. As 
mentioned by Coakley (ref. 15), the q-cj model was found to be remarkably insensitive to the 
initial field that is specified to start the calculations for q and w. In the present tests it has always 
been possible to start the model with a flat distribution of the turbulence quantities. For the 
viscous calculations, it was a good practice to keep CFL<10 together with 1/32 as a weight for 
the fourth-order artificial damping. The total absence of shocks allowed for eliminating the 
second-order artificial damping in all the calculations. Moreover, the aforementioned wall 
treatment of the fourth-order damping ensured inlet-outlet mass errors of the order of 0.5 percent. 

Concluding Remarks 

The scalar form of the approximate factorization coupled with a turbulence model proved to 
be suited for the solution of turbulent internal flows where diffusive terms play a dominant role. 
The introduction of an approximate treatment of the diffusive terms proved to increase 
convergence of the algorithm for internal flow configurations. Nevertheless, the implicit treatment 
of these terms must be tested for a wider range of geometries and Reynolds numbers. The 
influence of Mach number on convergence rate needs to be investigated, especially for the 
proposed approximate implicit treatment of diffusive terms based on a space derivative of the 
fluid density. New tests are currently being performed. 

For the Stanitz elbow geometry, the Baldwin-Lomax turbulence model proved to give results 
in acceptable agreement with the experiments, provided that the presence of multiple boundaries 
is properly accounted for. The low Reynolds number q-w two-equation model version investigated 
here proved to be suited for three-dimensional computations and gave a satisfactory description of 
the flow field with a manageable number of grid points and only a small increase in 
computational time with respect to the algebraic model. 
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RESIDUAL (L0G10) 


Table I. - Grids for the Stanitz elbow 


Grid 

number 

Points 

Expansion 

ratio 

1 

50x31x21 

1.2 

2 

50x41x31 

1.2 

3 

99x31x21 

1.2 

4 

50x31x21 

1.3 

5 

51x41x31 

1.3 


(a) no diffusive terms 


(b) equation 2.6 



(c) equation 2.7 
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Figure 1. - Convergence tests for implicit treatment of diffusive terms. 


(a) tridiagonal 


(b) pentadiagonal 


ITERATION N 


ITERATION .S 


Convergence tests for tri- or pentadiagonal solvers. 








log 10 DU 

- 5.0 - 4.5 - 4.0 - 3.5 - 3.0 - 2.5 - 2.0 


. . . i-m 


i-1 i-2 i-3 

i-4 

i-5 

i-6 

i-7 

00 

1 

•rt 

, * 

— * — 

— * — 

— * — 

— * — 

*-■ 

1 

+1 -4 +6 

-4 

i 

+1 




+1 -4 

1 

+6 

-4 

t 

+1 



+1 

-4 

1 

+6 

-4 

i 

+ 1 



+1 

-4 

1 

+6 

-4 

i 

+ 1 



+ 1 

-4 

1 

+ 6 

-4 

i 




+1 

-4 

1 

+6 





+ 1 

-4 






+1 


+1-3+3-10 0 0 0 


-1 i-2 i-3 

i-4 

i-5 

i-6 

| * * — 

* 



* 

1 

+1 -2 +1 




1 

+1 -2 

+ 1 
1 



+ 1 

1 

-2 

+ 1 
i 



+1 

1 

-2 

+ 1 
i 



+ 1 

1 

-2 


+1 -1 0 0 0 0 


* * 


(a) FOURTH ORDER 


. . . i— m 
* ★ 


(b) SECOND ORDER 


Figure 3. - Artificial dissipation treatment at boundaries. 
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Figure 5. - Sketch of experimental setup of S-duct. 
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Figure 7(b). - Velocity profiles on midspan of S-duct for laminar flow regime. 



Figure 7(c). - Secondary velocity at the bend exit of S-duct for laminar (low regime. 
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Figure 8(c). - Wall C p distribution of S-duct for turbulent flow regime. 



Figure 8(d). - Convergence history for turbulent flow regime. 
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Figure 11. - Fnviscid Ps distribution of Stanitz elbow. 



IMPLICIT 

IMPLICIT 


Figure 12. - Convergence history of Stanitz elbow. 




Figure 13. - Refined 50x41x31 point grid of Stanitz elbow. 
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Figure 14. - Ps distribution with multiple wall treatments in Baldwin-Lomax model. 
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Figure 21. - Secondary flow visualization of Stanitz elbow. 
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